clear
pr drop _all
set matsize 800
set more off

***********************
*  Likelihood Evaluator                                                           
***********************

program define splag_ll

args lnf mu rho sigma
tempvar A rSL
gen `A'= ones - `rho'*EIGS1
gen `rSL'=`rho'*SL1
qui replace `lnf'= ln(`A') + ln(normden($ML_y1-`rSL'-`mu', 0, `sigma'))
end

global nobs = 706

***********************
* Open Data For Weights
***********************
clear
use "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\govcon_w.dta", clear
***********************
mkmat var1-var$nobs, matrix(W)
matrix eigenvalues eig1 imaginaryv = W
matrix eig2 = eig1'
matrix ones=J($nobs,1,1)


**************************
* Open Data for Regression
**************************
drop _all
use "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\govcon.dta", clear
**************************
global Y gcgdp
global X gcgdp_1 lnimp_1 deind_1 imp_deind unemp imp_unemp c1-c19
*lnexp_1 rgdppc1 oldage1 leftc1 undens1 c1-c19 t2-t40
*unemp imp_unemp lnexp_1 rgdppc1 oldage1 leftc1 undens1 c1-c19 
*t2-t40
mkmat $Y, matrix(Y)
matrix SL = W*Y
svmat SL, n(SL)
svmat eig2, n(EIGS)
svmat ones, n(ones)

*************************
*Produce starting values 
*************************
qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

	ml model lf splag_ll (mu: $Y=$X) (rho:) (sigma:), technique(dfp) vce(opg)
	ml init OLSb
	ml init rho:_cons=0
	ml init sigma:_cons=`OLSsigma'
	*ml check
	*ml report
	ml max, difficult ltolerance(1e-8) tolerance(1e-8)

predict yhat
mkmat yhat, matrix(yhat)
matrix e2 = (Y-yhat)'*(Y-yhat) 
scalar e2 = e2[1,1]
qui summarize $Y
matrix ym = r(mean)*J($nobs,1,1)        
matrix ym2 = (Y-ym)'*(Y-ym)
scalar ym2 = ym2[1,1] 
scalar r2 = 1 - (e2/ym2)
scalar list r2


*******************Spatio-Temporal Response Path******************

***********************
* Open Data For Temportal Weights
***********************
drop _all
use "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\govcon_m2.dta", clear
***********************
mkmat var1-var$nobs, matrix(M)

scalar deind_1=60
scalar deind_2=75
matrix eye = I($nobs) 
matrix estcf1 = eye - _b[rho:_cons]*W - _b[mu:gcgdp_1]*M
matrix estcf_1 = (_b[mu:lnimp_1]+ _b[mu:imp_deind]*deind_1)*inv(estcf1)
matrix estcf_2 = (_b[mu:lnimp_1]+ _b[mu:imp_deind]*deind_2)*inv(estcf1)
matrix estcf_3 = _b[mu:imp_deind]*inv(estcf1)
matrix d_cf_rho1_1 = estcf_1*W*inv(estcf1)
matrix d_cf_rho1_2 = estcf_2*W*inv(estcf1)
matrix d_cf_rho1_3 = estcf_3*W*inv(estcf1)
matrix d_cf_rho_1 = d_cf_rho1_1[1..., 495]
**495 is Germany 1990**
matrix d_cf_rho_2 = d_cf_rho1_2[1..., 495]
matrix d_cf_rho_3 = d_cf_rho1_3[1..., 495]
matrix d_cf_phi1_1 = estcf_1*inv(estcf1)
matrix d_cf_phi1_2 = estcf_2*inv(estcf1)
matrix d_cf_phi1_3 = estcf_3*inv(estcf1)
matrix d_cf_phi_1 = d_cf_phi1_1[1..., 495]
matrix d_cf_phi_2 = d_cf_phi1_2[1..., 495]
matrix d_cf_phi_3 = d_cf_phi1_3[1..., 495]
matrix d_cf_beta1 = inv(estcf1)
matrix d_cf_beta = d_cf_beta1[1..., 495]
matrix d_cf_TEMP_1 = [d_cf_rho_1, d_cf_phi_1, d_cf_beta]
matrix d_cf_TEMP_2 = [d_cf_rho_2, d_cf_phi_2, d_cf_beta]
matrix d_cf_TEMP_3 = [d_cf_rho_3, d_cf_phi_3, d_cf_beta]
matrix d_cf_1 = d_cf_TEMP_1[495...,1...]
matrix d_cf_2 = d_cf_TEMP_2[495...,1...]
matrix d_cf_3 = d_cf_TEMP_3[495...,1...]
matrix V = e(V)
matrix VCVM_1 = [V[25,25], V[25,1], V[25,2]+deind_1*V[25,4] \ V[1,25], V[1,1], V[1,2]+deind_1*V[1,4] \ V[25,2]+deind_1*V[25,4], V[1,2]+deind_1*V[2,4], V[2,2]+(deind_1^2)*V[4,4]+2*V[2,4]]
matrix VCVM_2 = [V[25,25], V[25,1], V[25,2]+deind_2*V[25,4] \ V[1,25], V[1,1], V[1,2]+deind_2*V[1,4] \ V[25,2]+deind_2*V[25,4], V[1,2]+deind_2*V[2,4], V[2,2]+(deind_2^2)*V[4,4]+2*V[2,4]]
matrix VCVM_3 = [V[25,25], V[25,1], V[25,4] \ V[1,25], V[1,1], V[1,4] \ V[25,4], V[1,4], V[4,4]] 
matrix delta1_1 = d_cf_1*VCVM_1*d_cf_1'
matrix delta_1 = vecdiag(delta1_1)'
matrix delta1_2 = d_cf_2*VCVM_2*d_cf_2'
matrix delta_2 = vecdiag(delta1_2)'
matrix delta_3 = d_cf_3*VCVM_3*d_cf_3'
matrix delta_diff = vecdiag(delta_3)'
matrix cf_1 = estcf_1[495...,495]
matrix cf_2 = estcf_2[495...,495]  
matrix cf_3 = estcf_3[495...,495]
matrix cf_diff = cf_1 - cf_2
drop _all
svmat cf_1, n(cf_1)
svmat delta_1, n(delta_1)
svmat cf_2, n(cf_2)
svmat delta_2, n(delta_2)
svmat cf_diff, n(cf_diff)
svmat delta_diff, n(delta_diff)
outfile cf_11 delta_11 cf_21 delta_21 cf_diff1 delta_diff1 using "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\cf_response_path", wide replace

*************Long-Run Steady-State Spatio-Temporal Effects for deind=deind_2********************
* Open Data For NxN Weights
***********************
drop _all
use "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\govcon_w_LR.dta", clear
***********************
mkmat var1-var20, matrix(W_LR)
matrix eye_n=I(20) 
matrix estcf1_LR = eye_n - _b[rho:_cons]*W_LR - _b[mu:gcgdp_1]*eye_n
matrix estcf_LR = [_b[mu:lnimp_1] + _b[mu:imp_deind]*deind_2]*inv(estcf1_LR)
matrix d_cf_rho1_LR = estcf_LR*W_LR*inv(estcf1_LR)
matrix d_cf_rho_LR = d_cf_rho1_LR[1...,1] 
matrix d_cf_phi1_LR = estcf_LR*inv(estcf1_LR)
matrix d_cf_phi_LR = d_cf_phi1_LR[1...,1]
matrix d_cf_beta1_LR = inv(estcf1_LR)
matrix d_cf_beta_LR = d_cf_beta1_LR[1...,1]
matrix d_cf_TEMP_LR = [d_cf_rho_LR, d_cf_phi_LR, d_cf_beta_LR]
matrix d_cf_LR = d_cf_TEMP_LR
matrix delta1_LR = d_cf_LR*VCVM_2*d_cf_LR'
matrix delta_LR = vecdiag(delta1_LR)'
matrix cf_LR = estcf_LR
svmat cf_LR, n(cf_LR)
svmat delta_LR, n(delta_LR)

*************Short-Run Spatial Effects for deind=deind_2********************
matrix estcf1_SR = eye_n - _b[rho:_cons]*W_LR 
matrix estcf_SR = [_b[mu:lnimp_1] + _b[mu:imp_deind]*deind_2]*inv(estcf1_SR)
matrix d_cf_rho1_SR = estcf_SR*W_LR*inv(estcf1_SR)
matrix d_cf_rho_SR = d_cf_rho1_SR[1...,1] 
matrix d_cf_beta1_SR = inv(estcf1_SR)
matrix d_cf_beta_SR = d_cf_beta1_SR[1...,20]
matrix d_cf_TEMP_SR = [d_cf_rho_SR, d_cf_beta_SR]
matrix d_cf_SR = d_cf_TEMP_SR
matrix VCVM_4 = [V[25,25], V[25,2]+deind_2*V[25,4] \ V[25,2]+deind_2*V[25,4], V[2,2]+(deind_2^2)*V[4,4]+2*V[2,4]]
matrix delta1_SR = d_cf_SR*VCVM_4*d_cf_SR'
matrix delta_SR = vecdiag(delta1_SR)'
matrix cf_SR = estcf_SR
svmat cf_SR, n(cf_SR)
svmat delta_SR, n(delta_SR)
outfile cf_LR1 delta_LR1 cf_SR1 delta_SR1 using "C:\Documents and Settings\jchays\Desktop\FINAL_BOOK_STUFF\Chapter 2\cf_LR_SR", wide replace



*Alternative Likelihood: ln(`A') - .5*(ln(2*_pi)) - .5*ln((`sigma')^2) - ((.5/((`sigma')^2))*(($ML_y1-`rSL'-`mu')^2))
